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ABSTRACT 

We continue our deep optical imaging survey of the Virgo cluster using the CWRU Burrell Schmidt tele- 
scope by presenting B-band surface photometry of the core of the Virgo cluster in order to study the cluster's 
intracluster light (ICL). We find ICL features down to iib^ 29 mag arcsec"^, confirming the results of Mihos 
et al. (2005), who saw a vast web of low-surface brightness streams, arcs, plumes, and diffuse light in the Virgo 
cluster core using V-band imaging. By combining these two data sets, we are able to measure the optical colors 
of many of the cluster's low-surface brightness features. While much of our imaging area is contaminated by 
galactic cirrus, the cluster core near the cD galaxy, M87, is unobscured. We trace the color profile of M87 
out to over 2000", and find a blueing trend with radius, continuing out to the largest radii. Moreover, we have 
measured the colors of several ICL features which extend beyond M87's outermost reaches and find that they 
have similar colors to the M87's halo itself, B-V^ 0.8. The common colors of these features suggest that the 
extended outer envelopes of cD galaxies, such as M87, may be formed from similar streams, created by tidal 
interactions within the cluster, that have since dissolved into a smooth background in the cluster potential. 

Subject headings: galaxies: clusters: individual (Virgo) — galaxies: individual (M87) — galaxies: interactions 
— galaxies: photometry 



1. INTRODUCTION 

Massive galaxy clusters are known to contain a population 
of stars which reside outside of any of the cluster's galaxies, 
often referred to as intracluster light or ICL. ICL features typi- 
cally have extremely faint surface brightnesses of < 1% of the 
brightness of the night sky, making their study extremely dif- 
ficult. While the first indications of the existence of ICL came 
from observations by Zwicky (1951), only with the advent of 
modern CCD technologies have detailed studies of these stars 
been made possible (e.g., Uson et al. 1991; Vflchez-Gomez et 
al. 1994; Bernstein et al. 1995; Gregg & West 1998; Trentham 
& Mobasher 1998). 

The most straightforward method of detecting the ICL is 
through deep broadband imaging at optical wavelengths (e.g., 
Uson et al. 1991; Vflchez-Gomez et al. 1994; Trentham & 
Mobasher 1998; Feldmeier et al. 2002, 2004a; Mihos et al. 
2005, hereafter refeiTed to as M05; Gonzalez et al. 2005; 
Krick & Bernstein 2007). With such imaging, we can not 
only measure the luminosity of ICL the component, but we 
can study its spatial distribution and detect individual ICL fea- 
tures, such as streams, arcs, and plumes (e.g., Gregg & West 
1998; Trentham & Mobasher 1998; Calcaneo-Roldan et al. 
2000; White et al. 2003; M05; Krick et al. 2006; Yagi et al. 
2007) as well as any large-scale diffuse components (e.g., 
Gonzalez et al. 2000; Feldmeier et al. 2002, 2004a; Adami 
et al. 2005; Zibetti et al. 2005; Patel et al. 2006; Krick & Bern- 
stein 2007; Pierini et al. 2008; Da Rocha et al. 2008). ICL has 
also been detected using discrete stellar tracers (e.g., Fergu- 
son et al. 1998; Feldmeier et al. 1998; Durrell et al. 2002; Gal- 
Yam et al. 2003; Arnaboldi et al. 2004; Feldmeier et al. 2004b; 
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Gerhard et al. 2005; Aguerri et al. 2005; Neill et al. 2005; 
Maoz et al. 2005; Williams et al. 2007; Castro-Rodriguez 
et al. 2009; McGee & Balogh 2010), although these methods 
often lack the ability to detect individual ICL features. 

Intracluster light is thought to form primarily by the tidal 
stripping of stars as galaxies interact and merge during the 
hierarchical accretion history of the cluster, causing the frac- 
tion of the cluster's luminosity found in the ICL to increase as 
it evolves (e.g., Napolitano et al. 2003; Murante et al. 2004; 
Willman et al. 2004; Rudick et al. 2006; Monaco et al. 2006; 
Conroy et al. 2007; Murante et al. 2007; Purcell et al. 2007; 
Yan et al. 2009; Baria et al. 2009). Numerous mechanisms 
for generating the ICL have been proposed, including the in- 
fall of groups into the cluster potential (Willman et al. 2004; 
Rudick et al. 2006), high speed encounters within the cluster 
(Moore et al. 1996; Gnedin 2003), and galactic mergers dur- 
ing the buildup of the massive central galaxy (Murante et al. 
2007; Conroy et al. 2007). All of these processes are likely 
to be occurring simultaneously and each will create distinct 
observable signatures (Rudick et al. 2009). Thus, the forma- 
tion of ICL is intimately linked to the dynamical history of the 
cluster, and the observable features of the ICL should contain 
a great deal of information about the evolutionary processes 
which have shaped both the cluster and its constituent galax- 
ies. 

While ICL is generally thought of as stellar material found 
outside of any individual galaxy, in practice galaxies have no 
well-defined edge (Abadi et al. 2006). Thus, the distinction 
between the ICL and the outer luminosity profiles of cluster 
galaxies, which display similar surface brightnesses in broad- 
band imaging, is difficult and somewhat arbitrary (M05). In 
fact, simulations of the formation of clusters' most massive el- 
liptical galaxies have shown that their extended stellar profiles 
form through similar merger and tidal stripping mechanisms 
as the more diffuse ICL (Dubinski et al. 1998; Monaco et al. 
2006; Com-oy et al. 2007; Murante et al. 2007; Ruszkowski & 
Springel 2009). We therefore prefer to use intracluster light 
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as a qualitative description of cluster luminosity at low surface 
brightness, which may refer to individual tidal streams, ex- 
treme galactic outskirts, or any large-scale diffuse luminosity 
component, and which are all likely products of tidal stripping 
and disruption of galaxies during the dynamical evolution of 
the cluster. 

In addition to the quantity and morphology of the ICL, 
knowledge of the underlying stellar populations also provides 
a vital tool for understanding its formation. Giant elUptical 
galaxies are known to display radial color gradients, whereby 
the color index decreases, or becomes bluer, with increas- 
ing radius (e.g., Pettit 1954; de Vaucouleurs 1961; Carter & 
Dixon 1978; Strom & Strom 1978; Davis et al. 1985; Vader 
et al. 1988; Goudfrooij et al. 1994; Bernardi et al. 2003; 
Cantiello et al. 2005; Liu et al. 2005), primarily due to gradi- 
ents in the stellar metallicities (e.g., Spinrad et al. 1972; Strom 
et al. 1976; Baum et al. 1986; Carollo et al. 1993; Tamura 
et al. 2000; Sanchez-Blazquez et al. 2007; Rawle et al. 2010). 
Just as the strength of these gradients is an important clue in 
reconstructing the formation history of elliptical galaxies, the 
ages and metallicities of the intracluster stars will be highly 
dependent on the progenitor galaxies from which they were 
stripped. The results of Sommer-Larsen et al. (2005) show 
that the metallicity of the ICL is expected to be on average 
similar to that of the outer envelope of the cD galaxy, while 
Purcell et al. (2008) show that individual ICL streams formed 
by recent interactions should be more metal rich than the 
surrounding diffuse component. Murante et al. (2004) have 
shown that the diffuse ICL stellar population is expected to be 
older, on average, than the galactic stars. 

Observationally, it is very difficult to break the age- 
metallicity degeneracy and make precise determinations of ei- 
ther quantity. Although Williams et al. (2007) used HST ACS 
imaging of red giant branch stars in the Virgo cluster to do so, 
such observations are only possible in the most nearby sys- 
tems and span very small fields of view. Their work found that 
the ICL is composed of stars with a wide variety of ages and 
metallicities, but that the dominant component is old ( ^ 10 
Gyr) and moderately metal poor ([M/H] ^ - 1.0). A more 
common technique is to instead use broadband optical colors 
in order to compare the ICL stellar populations to those of 
galactic stars. Most studies have found results consistent with 
Williams et al. (2007), where the average ICL color is sim- 
ilar to that of the outer halo of the cluster's brightest galaxy 
(e.g., Zibetti et al. 2005; Pierini et al. 2008; da Rocha el al. 
2008), which due to radial color gradients is bluer than the 
galaxy's interior regions. However, a number of clusters have 
been observed to have an ICL component with significantly 
redder colors, more similar to the brightest galaxy's interior 
(Krick & Bemstein 2007; da Rocha & Mendes de Oliveira 
2005; Gonzalez et al. 2000). 

As part of our ongoing deep imaging survey of the Virgo 
cluster using Case Westem Reserve University's Burrell 
Schmidt telescope, this paper presents deep imaging of the 
Virgo cluster core in the B-band. By combining this data 
with the V-band results previously published in M05, we are 
able to measure the optical colors of Virgo's intracluster light. 
A detailed description of our data acquisition and reduction 
techniques is given in Section 2. Our B-band image is pre- 
sented in Section 3, while Section 4 combines the two images 
in order to measure the colors of ICL features. Section 5 in- 
cludes a summary of our results and a discussion of our inter- 
pretations. Finally, detailed error models for our photometric 
measurements can be found in the Appendix. 



2. OBSERVATIONS AND DATA REDUCTION 

2.1. Observations 

Our imaging data was collected on dark photometric nights 
on three separate observing runs spanning January through 
March of 2009, using the 0.6m CWRU Bun-ell Schmidt Tele- 
scope located at Kitt Peak National Observatory. The Burrell 
Schmidt's closed-tube optical design makes it ideal for wide- 
field, deep optical imaging. Using a 4096 x 4096 pixel CCD 
read out of 4 amplifiers simultaneously, the 1 .45" pixel scale 
results in a field of view which exceeds 2.5 square degrees. 
Further details of the CWRU Burrell Schmidt's optical sys- 
tem can be found in Nassau (1945) and Slater et al. (2009). 

The observations described in M05 were made using the 
Washington M filter, converted to Johnson V magnitudes. To 
measure the colors of the ICL features seen in those observa- 
tions, we could have selected to re-observe the area in either a 
bluer or redder filter. However, at redder wavelengths, emis- 
sion lines in the atmosphere cause the night sky to vary con- 
siderably in brightness on timescales of minutes (Feldmeier 
et al. 2002 and references therein), which would severely 
hamper our ability to accurately flat field and sky-subtract our 
data. While both our targets and the sky are fainter in the blue, 
requiring longer observing times, the stability of the night sky 
emission made the selection of a blue filter the natural choice. 

Because the Washington M filter is somewhat shifted blue- 
ward of the V filter, and to ensure adequate wavelength sepa- 
ration between our two filters, for these observations we used 
a custom-designed filter, referred to as B', with an effective 
wavelength « 200 shorter than the standard Johnson B fil- 
ter. All surface brightness measurements, however, have been 
transformed to the Johnson B system. 

Our observational and data reduction techniques are based 
on those described in Morrison et al. (1997), Feldmeier et al. 
(2002), M05, and Slater et al. (2009). Because of the diffuse, 
extremely faint nature of the ICL features we are measuring 
(0. 1 — 1 % of the brightness of the night sky), generating a flat- 
field which exceeds this precision over the entire field of view 
is of critical importance. We therefore constructed a night- 
sky flat from pointings of pre-selected sky fields, chosen to 
be offset from our target field by « 0.25 - 1 hours in right as- 
cension and less than « 4° in declination, while containing 
relatively low stellar density and no particularly bright stars. 
Our observing schedule was roughly evenly divided between 
imaging our target and the offset sky fields. A typical obser- 
vation cycle consisted of 2 sky images offset west, 4 images 
on the target field, and 2 more sky images offset east, repeated 
throughout the night. By minimizing the motion of the tele- 
scope during the cycle, this pattern reduced any hysteresis ef- 
fects caused by flexure of the optical system and ensured that 
the sky frames were taken as close in time and position as 
possible to the target images, under similar atmospheric and 
telescope conditions. When imaging the target field, individ- 
ual pointings were intentionally offset in a dither pattern from 
the field center by as much as 0.75°, or half the size of the 
field of view, in order to reduce systematic effects which may 
result from consistently imaging objects with the same area 
of the CCD. All sky and target field images were 1200 second 
integrations, yielding a sky brightness of ~ 750 ADU. 

2.1.1. Sources of Diffuse Light 

Although the scientific aim of our survey is to image the in- 
tracluster light in the Virgo galaxy cluster, we are, of course, 
sensitive to all sources of low-surface brightness diffuse light. 
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In practice, we detect two major sources of diffuse light: 
ICL in the Virgo cluster and light scattered by galactic cirrus 
within our own Galaxy (e.g., Sandage 1978; Guhathakurta & 
Tyson 1989; Witt et al. 2008). For the purposes of image 
processing and data reduction, all diffuse light sources are 
equivalent and we make no distinction between luminosity 
from galactic or extra-galactic sources. For scientific anal- 
yses, however, the source of the diffuse light can be critical 
to the interpretation of our results, and the impact of galactic 
cirrus is discussed in Section[J!2l 



2.2. Image Pre-Pmcessing, Photometric Solutions, and Data 
Quality Cuts 

Our image processing procedure began by subtracting a 
nightly bias frame, applying a CCD overscan region correc- 
tion, and removing amplifier crosstalk from all images, each 
in the usual manner (hereafter, these images are referred to 
as pre-processed). The images were then flattened using a 
preliminary flat field generated by a simple median combi- 
nation of all of our sky frames scaled by their mode^sing 
the imcombine task from the IRAF software package Fl Al- 
though the input sky images were not put through the rig- 
orous screening a nd pr ocessing required to create a final flat 
field (see Section 2.3 1, pixel values in this preliminary flat 
field, created from 111 individual images, differed from the 
final version by typically < 0.5%. An astrometric solution for 
each image was calculated using the Astrometry.net software 
package (Lang et al. 2010), and an airmass correction was ap- 
plied. 

Photometric zeropoints for all images were then calculated 
by comparing the stellar fluxes, as measured using IRAF's 
daophot routines, to those found in the Sloan Digital Sky Sur- 
vey DR7 catalog (Abazajian et al. 2009), transformed to John- 
son B magnitudes (Ivezic et al. 2007). From this analysis, we 
calculated a B-V color term of 0.05 mag to transform our 
data to the Johnson B system. For each night of observations 
we calculated and subtracted a nightly zeropoint, taken to be 
the mean of all zeropoints from that night. After subtracting 
this nightly zeropoint, all images were set to a common pho- 
tometric zeropoint. For a typical B-V color of 1 .0, 1 ADU 
corresponds to 29.6 mag arcsec"^. 

A data quality cut was made by excluding images from 
nights found to have highly variable photometric zeropoints, 
indicative of light cloud cover or other non-photometric con- 
ditions. We calculated the mean and scatter in zeropoint for 
all images and excluded any data from nights on which the 
photometric zeropoint varied by > 2a from the mean, or ap- 
proximately 0.02 mag. This quantitative measure of photo- 
metric atmospheric conditions agreed remarkably well with 
observer notes from the nightly observation logs. 

2.3. Flat Field 

Our flat field was constructed using the offset sky images 
described in Section |2.1[ The primary motivation for using 
this sky-flat method is that we expect the detector response to 
be wavelength-dependent. By using the night sky as our uni- 
form illumination source, we ensure that the flat field is mea- 
sured at the same effective wavelength as our target sources. 
Because the diffuse light sources we are measuring are < 1% 

^ IRAF is distributed by the National Optical Astronomy Observatory, 
which is operated by the Association of Universities for Research in Astron- 
omy (AURA) under cooperative agreement with the National Science Foun- 
dation. 



of the brightness of the night sky, we are dominated by sky 
photons even in our target frames. Moreover, we expect the 
ICL itself to be similar in color to the night sky (Taylor et al. 
2004; Sommer-Larsen et al. 2005), further motivating the use 
of the night sky as our uniform illumination source. 

The first step in creating the final flat field was to use 
IRAF's objmask task to find and mask all objects in each im- 
age, including stars, galaxies, satellite trails, and any other 
bright features. Because our goal is to measure a uniform il- 
lumination pattern as well as possible, we used an aggressive 
masking technique whereby we ran objmask twice, first on 
the original image and again after masking out any objects 
found in the first run. Once the object masks were applied, 
we median-binned the images in 128 x 128 pixel groups, and 
visually inspected each frame to check the image quality and 
search for and mask any remaining diffuse, non-sky features. 
All images had some features which were masked by hand, 
including the extended wings of bright stars, internal reflec- 
tions, and flares from stars near the edge of the field of view. 
A number of images were discarded due to the presence of 
strong galactic cirrus, a particularly bright star on the image, 
or other features which created a large-scale illumination pat- 
tern. After this quality cut we were left with 84 night sky 
images with which to create the final flat field. 

Due to variations in the atmospheric conditions of each ob- 
servation, such as airglow, transparency, airmass, etc., the 
mean flux level of the night sky varies from image to image, 
typically in the range of 650-850 ADU, or /iB=22.6-22.3 mag 
arcsec"^. However these effects not only cause the mean sky 
level to vary, they are also the cause of an illumination gra- 
dient over the field of view of each image, typically in the 
range of 5-10 ADU, or « 1% of the sky flux across the frame. 
In order to account for these atmospheric effects and create 
a flat field representing a uniform illumination pattern, we 
employed an iterative process to combine the night sky im- 
ages into our final flat field, following Morrison et al. (1997). 
We began by flattening the pre-processed sky images using 
the preliminary flat field and then applied both the object 
and hand masks. The sky gradients were removed by bin- 
ning and fitting a plane to each image, and then dividing the 
pre-processed images by this normalized sky plane. We then 
re-applied the object and hand masks, and median combined 
all the frames to make the new flat field. This process was 
then repeated, except that the pre-processed data were flat- 
tened with the new flat field. After only a few iterations, the 
resulting flat field converged. As a further check of the ro- 
bustness of our combination technique, we made several flat 
field images from a random selection of half of our sky im- 
ages; pixels varied in these half-split flats with a standard de- 
viation of sa 0.5%, consistent with the errors expected from 
photon statistics (see the Appendix for a discussion of photon 
statistics and other sources of photometric error in the final 
images). All pixels which were either under-illuminated by 
over 50% from the mean or which varied strongly in the half- 
split flats — these pixels were typically associated with small 
dust spots or very near the edge of the CCD, and represent 
< 0.5% of the pixels — were masked in the target frames. 

During our processing of the sky images into a flat field, we 
noticed a systematic difference between images taken during 
the third observing run and those taken during the first two 
runs. Frames from the third run showed a characteristic ra- 
dial dimming pattern whereby pixels farthest away from the 
image center were under-illuminated by 2%. We believe 
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Fig. 1. — Solid line: the extended wings of our model stellar profile for 
a mg = star. Dotted line: Our model PSF, excluding reflections from the 
dewar window surfaces. Dashed and dot-dash lines: the contributions to the 
stellar profile from the inner and outer dewar window surfaces, respectively, 
for a star with axisymmetric reflections. Unlike the data used in Slater et al. 
(2009) taken using the Washington-M filter, we see no measurable reflections 
from the B' filter surfaces. 

this effect to be the result of a thermal gradient, whereby the 
center of the CCD was cooled less efficiently, and was there- 
fore warmer, than the outskirts during this run. However, we 
find no indication that the flat field characteristics varied tem- 
porally within the third run. Fortunately, because the weather 
during the final run was significantly better than in the preced- 
ing runs, a disproportionate fraction of our data was taken dur- 
ing run three, giving us an ample number of exposures with 
which to construct a separate flat field. We thus generated two 
final flat fields, one for the first two runs using 47 sky images, 
and another for the third run from 37 sky frames. Our target 
images were then flattened by dividing by the final flat field 
corr esponding to the run in which they were taken. In Section 
|2.6| we show that the final mosaic image is not affected by any 
run-specific effects. 

2.4. Star Subtraction 

A key element of our data processing technique was to re- 
move scattered light from stellar sources in our target im- 
ages using the methods developed in Slater et al. (2009). As 
Slater et al. (2009) showed, the combination of extended stel- 
lar profiles plus scattered light from internal reflections within 
the optical system can cause a spurious diffuse signal which 
would contaminate the very faint ICL signal we are trying to 
measure. Moreover, the internal reflections are not axisym- 
metric about the star, and are dependent on the star's position 
within the field of view. Thus, in addition to taking both our 
sky and target images, during our observing runs we also took 
numerous images of bright stars, varying the position of the 
star within the field of view, in order to trace the stellar PSF 
and internal reflections. We acquired 13 900 s exposures of /3 
Gem {niB = 2.14), five 900 s exposures of a Boo (mg = 1.19), 
and four 1200 s exposures of a Gem (mg = 2.00), as well as 
several short (1-15 s) exposures of these stars. Using these 
images, we determined the amplitude and position of reflec- 
tions off the optical surfaces and measured the stellar PSF to 
a radius of over 0.5°. We show the extended wings of our 
model stellar profile in Figure [T] 



In order to subtract the stellar contribution from our tar- 
get images, for each image we created a map of the stellar 
light by reconstructing the profile of each star in the image, 
including the position-dependent internal reflections. Stellar 
brightnesses were measured from short (1-60 s) exposures of 
our target field taken during the observing runs and the pro- 
file of each star was scaled to match its flux. We mapped the 
profile of each star out to the radius at which it fell below 0.3 
ADU and subtracted these profiles from the individual target 
frames. 

2.5. Sky Subtraction 

In exactly the same manner as the sky frames discussed in 
each target image contains a sky flux level which 
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varies between images, and spatially across each image. How- 
ever, whereas fitting a plane to the sky signal was straightfor- 
ward for the sky images because the vast majority of their pix- 
els were of empty sky, the target field is permeated by large 
galaxies and extended diffuse light, meaning that there is es- 
sentially no true sky to measure. Because the sky and diffuse 
light signals vary by similar amplitudes over similar scales, it 
is extraordinarily difficult to disentangle the two and remove 
only the sky signal. 

We began the sky subtraction process by running two 
rounds of objmask, masking image artifacts such as star flares 
and reflections by hand, and spatially binning the image, in 
the same manner as we processed the sky images in Section 
|2.3| M05 attempted to use an iterative process to minimize the 
frame-to-frame scatter in the intensity of pixels at the same lo- 
cation after sky subtraction. However, subsequent testing has 
revealed that our implementation of this technique achieved 
results not significantly better than simply fitting a single sky 
plane to each target image (see the Appendix for details). We 
have therefore utilized this simpler technique, and the sky sig- 
nal was removed by fitting a plane to each target frame in the 
same manner as for the sky images. 

This inability to accurately fit the night-sky signal results 
in our largest source of systematic uncertainty across large 
scales, creating an uncertainty of 1-2 ADU across the major- 
ity of the image, discussed in detail in the Appendix. Ad- 
ditionally, because the entire target region is permeated with 
diffuse light, determining the sky zero level, or the flux level 
of "pure" sky, was also extremely difficult. We selected sev- 
eral regions of the image which looked by eye to be the dark- 
est, or least diffuse light-contaminated, and defined the sky 
level to be the mean pixel value in these regions. We empha- 
size that this is likely to be somewhat brighter than the true 
sky level as these regions probably contain some diffuse light, 
and that this could lead to a systematic over-subtraction of any 
large-scale diffuse light component. Fundamentally, these are 
issues which are inherent to any surface photometry measure- 
ments in which the objects fill the field of view, making a 
precise determination of the sky flux extremely difficult. 

2.6. Creating The Final Mosaic 

After sky subtraction, the final processing step was to sim- 
ply apply a photometric airmass correction to each target im- 
age. The fully-processed frames were then registered and all 
102 images were median combined into a mosaic using im- 
combine, using a 3-(j rejection to remove outlying pixel val- 
ues. All regions of the mosaic which did not contain a mini- 
mum of five exposures were masked. After creating the initial 
mosaic, a final sky plane was fit and removed, in the same 



5 



manner as described in Section 2.5 this plane had an ampli- 
tude of w 0.5 ADU across the width of the mosaic. 

As a test of the robustness of our sky subtraction and image 
mosaicking procedures, we created numerous mosaics from 
subsets of input images and compared these to the final mo- 
saic. The input frames were divided into subsets based on 
various parameters, including the run, airmass, hour angle, 
right ascension, declination, sky brightness, and local time of 
the observation. Only when we split the input images by right 
ascension (i.e. the position of the image within the mosaic), 
did we see any significant large-scale variations from the fi- 
nal mosaic, and these variations were of only w 0.5 ADU. 
This effect is most likely caused by specific diffuse light fea- 
tures which were present in certain areas of the mosaic which 
would have systematically influenced the fit of the sky plane. 

3. B-BAND MOSAIC IMAGE 

3.1. ICL Features 

Figure l2] shows our final mosaic, binned in 16x16 pixel 
groups (23.2" X 23.2"), with bright stars masked. Qualita- 
tively, the B-band image of the Virgo cluster core shown in 
Figure [2] appears almost exactly like the V-band image from 
M05, albeit with a significantly larger mosaic area due to the 
increased CCD size. In fact, every feature identified in M05 
appears in this image, including the distorted outer halo of 
M87, the two tidal streams extending northwest of M87, the 
doglegged plume to the north of NGC 4435/4438, the com- 
mon envelope around NGC 4413/IC 3363/IC 3349, and many 
other features as well (see M05 Figure 3 for a schematic find- 
ing chart of ICL features in this field). The identification of 
these features in a second imaging survey, using a different 
filter and with a modified instrumental set-up confirms the ro- 
bustness of our imaging techniques and the veracity of these 
extremely faint structures. 

3.2. Galactic Cirrus 

In addition to ICL in the Virgo cluster, another significant 
large-scale, low-surface brightness astrophysical feature seen 
in this image is galactic cirrus. Galactic cirrus is composed of 
dust clouds within the Galaxy which reflect galactic starlight 
(e.g., Sandage 1976; Witt et al. 2008). These cold dust clouds 
are most readily detected through their thermal emission in 
the far-infrared. Figure [3] shows the Schlegel et al. (1998) 
IRAS 100/im image of our mosaic field. The galactic cirrus 
features are most obviously seen in the southeast, northwest, 
and southwest corners of our Z?-band mosaic and match ex- 
tremely well to far-infrared emission. Fortuitously, the core 
of the Virgo cluster, which contains the majority of the ICL 
features we wish to measure, is relatively unaffected by galac- 
tic cirrus. However, the cluster core is surrounded by a ring 
of cirrus which obscures any ICL which may lie behind it. 
These galactic cirrus features were not as readily apparent in 
the V-band image of M05 since most of the cirrus lies outside 
of the smaller field of view of that data set. Because these 
cirrus features have no consistent optical color (Guhathakurta 
& Tyson 1989; Witt et al. 2008) and vary in intensity on very 
small scales, we do not attempt to model and subtract them 
from the data and must simply ignore any areas affected by 
these features for the purposes of measuring ICL. 

4. ICL B-V COLORS 

By combining our V-band imaging data from M05 with our 
current B-band dataset, we can measure the B-V colors of 



the diffuse light in the Virgo cluster core. However, because 
we have estimated the sky emission differently for our B-band 
data th an fo r the V-band data presented in M05 (described in 
Section [23| , for consistency we have re-reduced that data us- 
ing the new sky subtraction method. This has only a very 
minor effect on the resulting mosaic and all the features de- 
scribed in that work remain qualitatively similar; a more quan- 
titative discussion of the differences between these two sky 
subtraction techniques can be found in the Appendix. 

A major challenge in generating a large-scale color map of 
the cluster is setting the proper sky flux level for both images. 
Not only is the absolute image sky level difficult to determine 
precisely (see Section 2.5 1, but within each image we also 
have large-scale sky level uncertainties which we estimate to 
be of order 1-2 ADU (see the Appendix). Because the features 
we are measuring have very low flux levels, even small errors 
in the sky levels can translate into substantial differences in 
the final colors. For example, at /iy=28 mag arcsec"^, a sky 
level offset of just 1 ADU becomes a 0.4 mag. uncertainty 
in the B-V color At higher flux levels this effect is signif- 
icantly reduced as the sky level uncertainty becomes a very 
small fraction of the total flux; at /iv=26 mag arcsec"^, a 1 
ADU zero level offset translates to only a 0.07 mag. uncer- 
tainty in the B-V color Figure|4]shows our B-V color map in 
the region of overlap between the two mosaic images, which 
is limited primarily by the V-band areal coverage. The image 
is shown on the same pixel scale as Figure|2] with all pixels at 
surface brightness iiv> 27.0 mag arcsec"^r iib> 27.5 mag 
arcsec"^ masked. 

Each of the three giant elliptical galaxies — M87, M86, 
and M84 — shows a very red central core with a distinct ra- 
dial gradient trending bluer in the outskirts, in excellent agree- 
ment with previous studies (e.g.. Carter & Dixon 1978; Davis 
et al. 1985; Goudfrooij et al. 1994; Bernardi et al. 2003; Liu 
et al. 2005). In this image, however, each of these and several 
other galaxies shows an azimuthally varying color gradient in 
the innermost regions. This is an instrumental effect caused 
by very small (sub-arcsecond) offsets in the coordinate regis- 
trations of the B and V images, coupled with the very steep 
luminosity gradients present in the galactic centers. The two 
data sets were taken several years apart, with different CCDs, 
on a telescope which underwent substantial modification to its 
optical system in the intervening years. The resulting high- 
order distortions make it exceedingly difficult to register the 
entire field of view to the sub-arcsecond precision needed to 
precisely measure the rapidly varying inner regions of galax- 
ies, and our survey was never designed to accomplish this. 
Instead, our analyses are concentrated on the diffuse outer re- 
gions of galaxies which are not effected by these minute regis- 
tration offsets, due to the large angular sizes of these features. 

At very low surface brightnesses the sky level uncertainties 
discussed above make any determination of the color from 
such a large-scale map highly problematic. Indications of 
these effects can indeed be seen in the upper left corner of 
this image, particularly in the galaxies NGC 4473 and 4477, 
which appear particularly red. Due to the specific geometry 
of the images which compose our mosaic, this is the region 
in which our sky level uncertainties are largest. We believe 
our sky level uncertainties in this area to be up to 3-4 ADU, 
which could easily lead to uncertainties of several tenths of a 
magnitude in B - V color, even at the relatively high surface 
brightnesses shown. At surface brightnesses fainter than the 
/iv> 27.0 mag arcsec"^, fiB> 27.5 mag arcsec"^ limits shown 
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Fig. 2. — B-band mosaic image of the Virgo cluster core (nortli is up, east is left). The low-surface brightness features seen in this image are qualitatively 
extremely similar to those seen in the V-band imaging of M05 (see Figure 1). All regions containing fewer than five exposures have been masked. Also, stars 
have been masked out to the radius where their subtracted profile falls below 5 ADU. The pixels have been binned in 16x 16 groups (23.2"x23.2") in order to 
increase the signal-to-noise and display faint features. The centers of the three giant elliptical galaxies in the field, M87, M86, and M84 have been marked with 
X's, and the red lines indicate the /is=25 mag arcsec"- isophotes for a selection of the cluster's most luminous galaxies. 



in Figurep] these issues become even more pronounced, mak- 
ing such alarge scale map an unreHable source for determin- 
ing the color of ICL features. 

In order to better measure the colors of interesting ICL fea- 
tures, we have developed two methods with which to measure 
colors and estimate the photometric errors, depending on the 
features' angular scales. For large, degree-scale objects, such 
as the extended stellar envelope of M87, the dominant un- 
certainty comes from our sky level gradients, which we can 
estimate and include in our error model. Over the smaller, 
^ 10' scales of individual tidal streams, other systematic ef- 
fects dominate, and we have developed a local background 
subtraction technique which allows us to more robustly esti- 
mate and reduce systematic uncertainties. The sections below 
provide examples of these techniques and measure the colors 
of many of the image's most interesting low-surface bright- 
ness features, while the error budgets are described in the Ap- 
pendix. 

4.1. Extended Stellar Envelope ofM87 

Our analysis of the luminosities, shapes, other V-band pho- 
tometric properties of M87 and other Virgo giant ellipticals 
can be found in Janowiecki et al. (2010). Here, we focus 
solely on the B-V color properties of M87 . Because the outer 



regions of the other two giant elliptical galaxies in the field, 
M84 and M86, overlap one another in projection, we have not 
attempted similar analyses for these galaxies. 

We have measured M87's color in elliptical annular bins, 
defined by the galaxy's V-band isophotes as measured by 
Janowiecki et al. (2010). As in that work, before measur- 
ing the galaxy luminosity, small-scale sources were masked 
using objmask. To effectively detect these sources, we sub- 
tracted a ring-median smoothed image from the original, thus 
removing the large-scale galaxy luminosity distribution and 
leaving only the small-scale sources (see Janowiecki et al. 
2010 for details). All sources detected in either image were 
masked. Figure [5] shows our resulting color profile. The open 
squares show the mean color within each annulus. The pho- 
tometric errors in the color measurements are entirely dom- 
inated by the systematic uncertainty in the sky level calibra- 
tion, which we estimate to be ±1 ADU on these scales (see the 
Appendix); the error bars in Figure |5]reflect this ±1 ADU un- 
certainty. Consistent with previous studies, our results show a 
bluing trend with radius. The dotted line in Figure [5] shows a 
fit of the color profile inside of 1000 arcsec, which has a slope 
of -0.11 in A(Z?-y)/Alog(5MA). Our measured color profile 
of M87's outer regions is in excellent agreement with the re- 
sults of Liu et al. (2005), and is consistent with a continuation 



Fig. 3. — IRAS 100/xm map of the Virgo cluster core from Schlegel et al. (1998), on the same pixel scale as Figure|2]with the B-band exposure mask applied. 
As in Figure |2] the centers of the three giant elliptical galaxies in the field. M87, M86, and M84 have been markeo with X's. and the red lines indicate the 
/is=25 mag arcsec"^ isophotes for a selection of the cluster's most luminous galaxies. Diffuse features are indicative of galactic cirrus, including the large knot 
in the northwest, the bright emission in the southeastern comer, and the long streak running nearly parallel to the southwestern border of the mosaic. We find 
con'esponding optical emission for all of these features, and must excl ude the se regions from our analysis of the ICL. Fortuitously, the region to the northwest of 
M87 containing the three prominent tidal features discussed in Section [4.2. 1 [ contains very little galactic cirrus. 

4.2. Colors of Streams 



of the color profile of the galaxy's inner regions measured by 
Zeilinger et al. (1993). 

To measure the azimuthal variance in the color, we have di- 
vided each elliptical annulus into eight equal angle slices, with 
the color results plotted as solid circles in FigurelS] At all radii 
beyond 400", the measured colors in all the angular slices lie 
within our ±1 ADU systematic error. Inside of this radius, 
the azimuthal variance is dominated by the coordinate regis- 
tration offsets discussed above. A closer examination of the 
measured colors in the outermost, lowest-surface brightness 
regions of the galaxy reveals that the color varies smoothly 
around the galaxy, indicating that the azimuthal variance is 
caused by a systematic gradient, such as that produced by sky 
level uncertainties, as opposed to random fluctuations in the 
photometry. Complicating these measurements even further 
is the presence of galactic cirrus which overlaps with M87's 
outermost reaches to the south and east. As described in Sec- 
tion 3.2 it is extremely difficult to separate the two compo- 



nents, and any photometric measurements of this area may, 
in fact, be dominated by the cirrus over M87's stellar halo. 
Removing the data from the most cirrus -contaminated areas, 
however, does not significantly improve the color scatter of 
the angular slices in Figure |5] 



While our ability to accurately measure colors is limited 
over large scales by gradients in the photometric sky levels 
and extended diffuse light sources, over shorter scales these 
uncertainties are much smaller Thus, for smaller scale fea- 
tures, such as individual ICL streams which span ^ 10', we 
can use regions adjacent to the features to define the local 
background flux and more precisely measure luminosities and 
colors. We have used a differential photometry technique 
whereby we define object and neighboring background re- 
gions, and measure the object luminosity by subtracting the 
flux in these local background fields. Discrete objects in the 
fields, such as bright galaxies, are masked so that we are mea- 
suring only the diffuse light of the features themselves, and 
we vary a number of parameters, including the precise place- 
ment of region boundaries, in order to estimate the measure- 
ment uncertainties (these processes are described in detail in 
the Appendix). In essence, this technique is a variation on 
familiar aperture photometry methods used for point sources, 
where the background flux is measured from an annulus just 
beyond a circular aperture. 

The background regions contain three sources of flux, each 
of which is also present in the object region and should be re- 
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Fig. 4. — B-V color map of the Virgo core in tlie regions of overlap between our two mosaic images, on the same pixel scale as Figure |2] As in Figure|2] 
the centers of the three giant elliptical galaxies in the field, M87, M86, and M84, have been marked with X's, and the red lines indicate the /is=25 mag arcsec 
isophotes for a selection of the cluster's most luminous galaxies. Bright stars and all pixels with ^v> 27.0 mag arcsec"- or ^b> 27.5 mag ai'csec"^ have been 
masked. As discussed in detail in the text, features at lower surface brightness are subject to large systematic errors which cause measurements of their colors 
to be unreliable over large scales. The colors in the central regions of galaxies have an asymmetric profile due to a small ( ^ 1") coordinate registration offset 
between the B and V images. 



moved to isolate the stream's luminosity. First, as discussed 
previously, the background regions should share any system- 
atic offset in the sky zero level with the object. Second, in ad- 
dition to individual tidal streams, the cluster ICL likely con- 
tains a large-scale diffuse component which permeates both 
the object and background regions. Finally, there are numer- 
ous foreground and background sources in all the regions; al- 
though we attempt to mask these objects (described in the Ap- 
pendix), because faint sources should be distributed approxi- 
mately homogeneously on these scales (Brainerd et al. 1995; 
Villumsen et al. 1997; Connolly et al. 2002; Coil et al. 2004; 
Morganson & Blandford 2009), their contribution should be 
equivalent in all regions. We discuss a detailed error budget 
for this photometric technique in the Appendix. 

4.2. 1 . M87's Northern Stellar Envelope and Tidal Streams 

A number of the most interesting low-surface brightness 
features in our mosaic occur in the northernmost regions of 
M87's extended stellar halo. In particular. Figure [6]highlights 
three ICL features that we have identified for more detailed 
study: the large tidal stream extending north-west from M87 
toward NGC 4461/4458 (stream A from M05); the smaller 
stream just to the west which is also emanating from M87 



TABLE 1 

Photometry of Low-Surface Brightness 
Features 



Region 


B-V 




(see Figures|6j&|^ 


3 


mag arcsec"^ 


A 


0.75-1.05 


28.6 


B 


0.8-1.2 


29.2 


C 


0.7-1.0 


28.7 


El 


0.8-0.9 


28.2 


E2 


0.75-0.85 


27.7 


E3 


0.75-0.85 


27.6 


Dl 


0.65-0.75 


27.2 


D2 


0.5-0.7 


28.4 



to the northwest (stream B from M05); and the broad plume 
to the east of these streams, almost directly north of M87's 
core (unlabeled in M05, labeled C in Figure |6]l. Addition- 
ally, we have defined three neighboring regions (labeled El, 
E2, and E3) which are more clearly part of M87's outer stel- 
lar envelope. In addition to being discussed below, all color 
measurements for these features can be found in Table [T] 
Precisely measuring the colors of these features is quite dif- 
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Fig. 5. — Color profile of M87's extended stellar halo. Squares are the 
color using all pixels in the bin, with ±1 ADU error bars. Circles show the 
azimuthal color variance by measuring the color in eight equal angle slices at 
each radius (see text for details), offset slightly right for clarity. The dotted 
line is a fit to the color profile inside of 1000", and has a slope of A(B — 
V)/A(log(SMA)) = -0.11. 

ficult, and the uncertainties are dominated by numerous sys- 
tematic effects which are detailed in the Appendix. For this 
reason, we prefer to quote our measured colors as ranges, 
rather than as central values with error estimates, which im- 
ply a well-characterized Gaussian error model. For features 
A, B, and C, we measure B-V colors of 0.75- 1.05, 0.8-1.2, 
and 0.7- 1.0, respectively, with mean surface brightnesses of 
28.6 mag arcsec"^, 29.2 mag arcsec"^, /iv~ 28.7 
mag arcsec"^, respectively. If these streams were displayed 
on the raw color map of Figure |4] however, they would each 
show a B-V color of ss 0.5. This difference in color high- 
lights the difficulty in interpreting the colors of very low sur- 
face brightness features from the raw B-V color map and 
the importance of using a local background subtraction when 
performing photometric measurements. 

Regions El, E2, and E3 in M87's stellar envelope display 
very similar colors to the tidal features. For El we measure 
B-V~ 0.8-0.9 with mean surface brightness of /iysa 28.2, 
while both E2 and E3 have B-Vw 0.75-0.85, with sur- 
face brightnesses of /lysa 27.7 mag arcsec~^and jiyK, 27.6 
mag arcsec"^, respectively. The measured color ranges are 
tighter in these regions due to their higher surface bright- 
nesses, meaning that systematic effects have a comparatively 
smaller effect (see the Appendix for details). Both E2 and 
E3 lie approximately along the 2000" semi-major axis ellipse 
around M87, and their colors measured using this differential 
method agree very well with those of the radial profile, shown 
in Figure |5] at this radius. 

Thus, within the measurement errors, we find that the op- 
tical colors of the tidal features are consistent with those of 
M87's outer stellar envelope. This suggests that these various 
features consist of similar stellar populations and may have a 
common origin. In fact, the colors we measure in these fea- 
tures are also consistent with the optical colors measured for 
Virgo's dwarf elliptical galaxy population of B-Vw 0.8 (van 
Zee et al. 2004). The implications of these findings for the 
evolution of galaxies within the Virgo cluster are discussed in 
Section|5] 



4.2.2. The Plume ofNGC 4435/4438 

Figure |7] shows another prominent low-surface brightness 
feature in the Virgo cluster core, the dog-legged plume just to 
the north of the interacting pair NGC 4435/4438 (feature D 
from M05). For this plume, we have divided the feature into 
two sections, the relatively high surface brightness vertical 
plume, and the lower surface brightness diffuse emission to 
the west, labeled Dl and D2, respectively, in Figure |7] These 
features are the bluest in the image for which we have made 
reliable measurements. We find that Dl has a B-V color of 
« 0.65-0.75 with mean surface brightness /iv« 27.1 mag 
arcsec"^ while D2 has B-V~ 0.5-0.7 and fiy^ 28.4 mag 
arcsec"^. 

A wide range of multi-wavelength observations have re- 
cently been made of this feature (Cortese et al. 2010; Krick 
et al. 2010), and its origin — whether a tidal feature from the 
interaction of NGC 4435/4438 or simply foreground galac- 
tic cirrus — is currently in doubt. While its visual mor- 
phology is immediately reminiscent of an interaction-induced 
tidal plume, especially given that NGC 4438 displays obvi- 
ous signs of tidal disturbance, Cortese et al. (2010) argue that 
the feature's UV-IR color and extremely narrow CO and HI 
velocity widths are more consistent with it being a galactic 
cirrus dust cloud. If it is truly an extragalactic feature, its blue 
colors may hint at recent tidally induced star formation, and 
HST imaging should be able to resolve a young stellar pop- 
ulation. If it proves to be a galactic cirrus feature, however, 
it provides an excellent, if sobering, example of the difficulty 
of discerning an object's origin by morphology alone and of 
the problems that galactic cirrus poses for deep observations 
of extragalactic objects. 

4.2.3. Limits of the Differential Photometry Technique 

While our differential photometry technique allows us to 
increase the precision of our color measurements by more 
robustly estimating the photometric uncertainties, the proce- 
dure's applicability is often limited by features' geometry. 
Most importantly, any features on which we use this tech- 
nique must have obvious background regions immediately 
surrounding them, preferably on multiple sides. As detailed 
in the Appendix, the choice of background region is often the 
dominant source of systematic uncertainty in a feature's mea- 
sured color This geometric criterion rules out many of the 
most interesting ICL features seen in our image as candidates 
for the differential photometry method. For instance, a re- 
gion of our mosaic displaying a number of interesting low- 
surface brightness tidal features is the area surrounding M86 
and M84, especially the many interacting galaxies to their 
south. However, because there are numerous overlapping ICL 
features throughout the region, including both individual tidal 
streams as well as larger diffuse components, it is extraordi- 
narily difficult to properly define regions in which to measure 
the background flux for any particular feature, thus limiting 
our ability to perform accurate photometry of these features. 

4.3. Fundamental Photometric Limits 

The extremely large color scatter in the azimuthally sliced 
bins of the outermost regions of M87's extended stellar halo 
seen in Figure [5] suggests that we are unable to measure the 
colors of diffuse light at these very low surface brightnesses 
with the precision necessary to meaningfully constrain the 
stellar populations present. We estimate the surface bright- 
ness limit to which we can adequately measure colors to be 
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Fig. 6. — The area just to the north of M87, containing several interesting tidal features for which we have measured B-V colors. The regions outlined in green 
- labeled A, B, and C, respectively — are the tidal features themselves. The blue regions — labeled El, E2, and E3 — are regions within M87's stellar envelope. 
The red regions are examples of background regions used to measure the local sky brightness. Discrete objects in the field, such as bright galaxies, are masked 
during the photometric analysis. Details of the photometric techniques used on these regions are given in the text. 

27 mag arcsec"^ in large, degree-scale features such as 
the giant elliptical galaxies. While our loca l background sub- 
traction technique detailed in Section |4!2] is able to substan- 
tially reduce the photometric errors for certain features with 
smaller angular sizes ( ^ 10 arcmin), allowing us to push our 
measurements to lower surface brightnesses, it is not univer- 
sally applicable to all diffuse light features in the image. Be- 
cause at the lowest surface brightnesses we are limited by sys- 
tematic uncertainties in the brightness of both the night sky 
and astrophysical sources of diffuse light, and not from ran- 
dom statistical noise or known features of the optical design, 
we may be approaching the fundamental limit of precision 
available from ground-based, wide-field surface photometry. 

While the Virgo cluster is an appealing target in which to 
study ICL in part because it is the nearest massive galaxy 
cluster and thus ICL features can be seen over large angular 
scales, its large angular size actually acts as a hindrance to our 
ability to perform accurate surfac e ph otometry at the faintest 
levels. As discussed in Section |2.5| even with the Burrell 
Schmidt's extremely wide-field imaging capability, the dif- 
fuse light features of the Virgo cluster essentially fill the en- 
tire field of view, making it extremely difficult to measure and 
subtract the brightness of the night sky. Thus, we may be able 
to achieve even greater imaging depth and reduce large-scale 
uncertainties for targets with smaller angular sizes, such as 
individual galaxies in the local universe or clusters at higher 
redshift, where the object does not fill the field of view, and 
the sky brightness, including large-scale sky gradients, can be 




Fig. 7. — The plume to the north of the interacting galaxy pair 
NGC4435/4438. Similar to Figure |6] the green areas — labeled Dl and 
D2 — enclose the plume itself, while the red regions are used as the local 
background regions. 
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more precisely measured. 

5 . SUMMARY AND DISCUSSION 

In this paper, we have presented our B-band image of the 
Virgo cluster core, taken as part of our ongoing deep imaging 
survey of the Virgo cluster All aspects of the survey, from the 
optical and mechanical design of the telescope system, to the 
acquisition, analysis, and reduction of the data have been op- 
timized to detect diffuse light at low surface brightness. Our 
final photometric uncertainties are dominated by systematic 
errors resulting from sky subtraction and astrophysical back- 
grounds and not simple photon statistics. 

Our Z?-band imaging confirms the results of M05, which 
found a vast web of diffuse light features in the core of the 
Virgo cluster These features include large-scale diffuse com- 
ponents, as well as a number of discrete features such as 
streamers, arcs, tails, and plumes which provide a wealth of 
information on the dynamical history of the cluster and its 
constituent galaxies. Just outside the cluster core, however, 
we have detected a large number of galactic cirrus features 
which prevent us from studying diffuse extragalactic sources 
which lie behind them. 

By combining these B-band results with the V-band imag- 
ing from M05, we have been able to measure the colors of 
a number of low surface brightness features in Virgo's core. 
We have measured the radial color profile of M87 out to very 
large radius, demonstrating that the color gradients seen in 
the inner regions of the galaxy extend to its outer stellar en- 
velope. Furthermore, we have measured the colors of several 
tidal features which extend from M87's stellar envelope and 
find that the two populations have similar optical colors within 
the measurement uncertainties. These results are consistent 
with the majority of other measurements of ICL color, which 
suggest that the ICL colors should be similar to those of the 
outskirts of the cluster's brightest galaxies (e.g., Zibetti et al. 
2005; Sommer-Larsen et al. 2005; Krick & Bernstein 2007). 

These results are consistent with the hypothesis that the 
outer envelopes of cD galaxies like M87 may be predomi- 
nantly built-up by the tidal stripping and disruption of smaller 
satellite galaxies, the same mechanism which likely generates 
the tidal streams and other ICL features (Rudick et al. 2009). 
In this scenario, the cD envelope is simply composed of tidal 



streams of ICL which have been mixed in the cluster potential, 
and dissolved to form a smooth distribution. While optical 
colors alone cannot definitively prove that any two stellar pop- 
ulations are identical due to the well-known age-metallicity 
degeneracy, the similar colors of these two populations is sug- 
gestive of a common origin. 

The optical colors of both the tidal features and M87's outer 
envelope are also consistent with the colors of the Virgo clus- 
ter's dwarf elliptical galaxy population (van Zee et al. 2004). 
Thus, the disruption of these galaxies in the cluster potential 
may provide a ready source for the intracluster stars which 
make up the low-surface brightness outer envelope of M87 
and the tidal streams. In fact, the thinness of two of the 



tidal streams studied in Section 4.2.1 suggests that they orig- 
inate from low-mass galaxies with small velocity dispersions 
(M05), such as dwarf ellipticals. Furthermore, the results of 
Williams et al. (2007), who measured the age and metallicity 
of Virgo's intracluster stars at much larger radius from any of 
the cluster's large galaxies to be ^10 Gyr and [M/H] ^ -1.0, 
are consistent with the same B-V^i 0.8 color that we mea- 
sure for the streams and stellar envelope immediately beyond 
M87 (Bmzual & Chai'lot 2003). This color similarity fur- 
ther bolsters the case that M87's extended stellar envelope, 
the tidal streams from disrupted galaxies, and the truly intra- 
cluster stars all form a single population, formed through sim- 
ilar mechanisms of tidal stripping during the cluster's hierar- 
chical assembly. 
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APPENDIX 
PHOTOMETRIC ERRORS 

There are a number of both random and systematic effects which contribute to the uncertainty in our photometric measurements, 
many of which were mentioned in the main body of the paper. Here, we discuss how these various effects interact to yield our 
final measurement errors. In general, we find that random errors from photon noise make an insignificant contribution to the 
error budget, and that the systematic effects which dominate are highly dependent on the angular scales being probed. Below, we 
calculate the errors from photon statistics to show that the y are negligible, and then describe the error sources which predominate 
in the angular scale regimes discussed in Sections [4. 1 |and[4!2] 



Small-scale Random Errors 

The pixel-to-pixel measurement error based simply on random photon statistics is a straightforward calculation. For a single 
image with a typical sky value of 750 ADU (because the ICL features we are measuring are < 1% of the sky brightness, we are 
dominated by sky photons even in our target frames), the photon noise in ADU is given by ^/Cs/g, where Cj is the sky flux in 
ADU and g is the gain in e" ADU"'. For our gain of 2.0e~ ADU"' this comes to 19.4 ADU pix"', or 2.6%. When we combine 
images, however, this error is scaled by a factor of 1.22/V5v, where is the number of images combined. Our flat field, which 
contains a minimum of 37 images, thus has an error of 0.52%, corresponding to 3.9 ADU pix"'. In the final mosaic, we further 
reduce the errors by combining a minimum of five images, for a final error of 0.28%, or 2.1 ADU pix"'. This is, however, an 
absolute maximum for our random photometric error, and for all measurements made in this paper it is significantly lower. 

In the calculation above we assumed the minimum number of five target exposures, whereas 30 exposures are more typical 
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Fig. 8. — The difference between the current image using a new sky subtraction technique and the 2005 image, in ADU. The butterfly pattern seen near the 
centers of bright galaxies is caused by a slight differences in the astrometric calibration of the two images (see Section|4j. 

across the image, which reduces the error by another factor of 2.6, or to < 1 ADU pix"'. Even greater reductions in the random 
errors, however, are achieved due to the fact that the diffuse nature of ICL features means that they are spread over a very large 
number of pixels. Any feature which covers even 100 pixels will have its random photometric error reduced by an order of 
magnitude. Because the final mosaic shown in Figure|2]has been binned in 16x16 pixel groups, its random per pixel photometric 
error from photon statistics falls to well below to 0.1 ADU, or < 0.01%. 

Large-scale ( ^ 0.5°) Systematic Sky Level Uncertainties 

As mentioned in Section [23] our ability to properly subtract the flux of the night sky from our object images leads to our 
greatest source of large-scale systematic uncertainty. In that section, we described the two different techniques that we have used 
to measure the sky signal. Figure IS] shows the difference between our V-band images, sky subtracted using the two methods. 
We expect the large-scale sky leveluncertainties in the V-band mosaic to be more severe than in the B-band, since the V-band 
data was taken using a CCD with only half of the current field of view, meaning that there were even fewer pixels available from 
which to estimate the sky level gradients. In general, the sky level offset varies smoothly across the field. The largest differences 
are found in the extreme northeast and southwest corners, where the two images differ by up to 4 ADU, but values of « 1-2 
ADU are much more common. While in certain specific cases, we see that one of the techniques yields unphysical results in 
certain areas of the image (e.g., the extremely red colors of NGC 4473 and 4477 seen in the northeast corner of Figure |4]l, we 
cannot constrain t he sk y level more generally. Thus, when measuring large-scale features, such as the extended stellar envelope 
of M87 in Section [4~T| we can only roughly estimate a systematic sky level uncertainty, using the variations seen in Figure[8]as a 
guide, and recognizing that this is not the only source of systematic uncertainty. 

Intermediate-scale ( ^ 10') Differential Photometry Uncertainties 

In Section |4.2| we described the differential photometry technique which we have developed to measure the colors of low- 
surface brightness features on intermediate scales, such as individual tidal features. Using this method, we can more robustly 
measure the colors of diffuse light features, resulting in much smaller uncertainties. We have identified four primary contributions 
to the uncertainty using this technique, each of which behaves in systematic, highly non-Gaussian ways. While the precise 
magnitude of these effects varies depending on the feature being measured, we use stream A, shown in Figure |6j as an example 
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to illustrate the relative importance of each. For this feature, we find a color of B — V~ 0.75- 1.05, which might be re-written as 
fi-y=0.9±0.15. 

In order to measure a feature's color, we need to determine the per pixel flux level of both the object and background regions. 
These measurements are limited by tiie uniformity of flux across the regions from astrophysical sources. While in the idealized 
case of a Gaussian distribution the precision to which we could determine the mean value is given by cr/ \/N, which for our data 
comes to sa 0.01 ADU, the non-Gaussianity of the pixel -values results in a much larger error We have estimated this uncertainty 
by using various statistical measures of the per pixel flux, including the assigning the per pixel flux to be the mean or median 
value. In general, we find that our results are relatively insensitive to the statistical measurement used, and for our example stream 
tills results in a B-V error of ± w 0.05 mag. 

There are two distinct effects contributing to tiie photometric uncertainty related to the placement of the object and backgroimd 
regions, respectively. The region boundaries were determined by eye, and are designed to qualitatively select the stream itself 
and the surrounding background regions where there are no obvious ICL features. In general, defining the object regions was 
relatively straightforward, and we find that almost any reasonable selection for the regions' boundaries have a comparatively 
small effect on the measurement uncertainty, which we estimate to be ± 0.05 mag in the B-V color for our example stream. 
The uncertainty resulting from selecting various background regions, however, has a profound effect on the measured colors, 
and is often the single largest source of error. Each background region has a slightly different flux level due to any small-scale 
sky level gradients or very faint diffuse light features present. Because the objects we are measuring have such low surface 
brightnesses, even small differences in the background flux levels result in large errors in the colors. For our example feature we 
estimate the uncertainty due to background region selection to be ± « 0.15 mag inB-V. This uncertainty from the selection of 
background regions stems primarily from using regions on different areas of the sky, whereas adjusting the precise borders of a 
backgroimd region has a comparatively smaU effect. 

The final component of the color uncertainty comes from our treatment of foreground/background sources which, of course, 
permeate the image. Ideally, any method we use to mask these sources should result in a negligible change in the color since 
it should affect all regions equally, and faint sources are expected to be unclustered on these scales (e.g., Brainerd et al. 1995). 
We have used two methods to eliminate the luminosity from bright sources. In the first, we ran objmask to mask all objects in 
the image; our final mask combined all objects detected in either the V or B images. Alternatively, we used only a cr-clipping 
algorithm to remove high-intensity pixels. We notice a systematic offset in the colors using these two masking procedures where 
the objmask technique produces redder B — V colors by « 0.1 mag. Unfortunately, using objmask to mask sources becomes 
unfeasible at flux levels much above 5 ADU (~ 27.5 mag arcsec"^), as the features we are measuring are themselves masked. 

Each feature's final quoted color range takes into account all of these effects by making numerous measurements while sys- 
tematically varying the input parameters. While we can estimate the uncertainties caused by each of these four effects, they are 
all systematic and interdependent sources of error, and thus not subject to the usual methods of error propagation. We quote the 
final measured colors as ranges to emphasize the non-Gaussianity of the errors, especially given the fact that our results are often 
not consistent with a central measured value. 

Although this paper is focused almost exclusively on the color of ICL features, and not their total luminosity or mean surface 
brightness, this differential photometry method is also suitable for measuring those quantities, and mean surface brightness is 
given for each feature in the text. One interesting detail that we find when we calculate these quantities is that the uncertainty in 
the color measurement is smaller than the uncertainty in mean surface brightness in either band; i.e., we can measure the ratio 
of the fluxes in the two bands more precisely than we can measure the absolute flux in either band. This peculiar effect arises 
from the fact that the dominant source of error in this photometric technique is the choice of background region. Because the 
variation between background regions is caused predominantly by changes in the flux level of astrophysical backgrounds, such as 
large-scale ICL components, the flux in these background regions is correlated in the two bands, and the mean surface brightness 
measurements are not independent of one another. This odd behavior of the photometric uncertainties once again illustrates 
that because our measurements are truly limited by astrophysical backgroimds, the resulting errors are highly systematic and 
non-Gaussian. 
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